﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.FFTPack
{
    /// <summary>
    /// Cosine FFT transform with odd wave numbers. 
    /// </summary>
    [Serializable]
    public class RealDoubleFFTEvenOdd : RealDoubleFFTMixed
    {
        /// <summary>
        /// This field holds the norm factor.
        /// </summary>
        private double normFactor;

        /// <summary>
        /// This field holds the wave table.
        /// </summary>
        protected double[] wavetable;

        /// <summary>
        /// This field holds the dimension.
        /// </summary>
        protected int ndim;

        /// <summary>
        /// Initializes a new instance of the <see cref="RealDoubleFFTEvenOdd"/> class.
        /// </summary>
        /// <param name="n">The size of a real data sequence. When n is a multiplication of 
        /// small numbers(4, 2, 3, 5), this FFT transform is very efficient.</param>
        public RealDoubleFFTEvenOdd(int n)
        {
            ndim = n;
            normFactor = 4 * n;
            if (wavetable == null || wavetable.Length != (3 * ndim + 15))
            {
                wavetable = new double[3 * ndim + 15];
            }
            Cosqi(ndim, wavetable);
        }

        /// <summary>
        /// This field can be used to normalize this FFT transform. This is because a call of forward 
        /// transform followed by a call of backward transform will multiply the input sequence by the 
        /// norm factor.
        /// </summary>
        /// <value>The norm factor to use for the transform.</value>
        public double NormFactor
        {
            get { return normFactor; }
            set { normFactor = value; }
        }

        /// <summary>
        /// Forward FFT transform of quarter wave data. It computes the coeffients in cosine series 
        /// representation with only odd wave numbers.
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x 
        /// contains the transform coeffients.</param>
        public virtual void ForwardTransform(double[] x)
        {
            Cosqf(ndim, x, wavetable);
        }

        /// <summary>
        /// Backward FFT transform of quarter wave data. It is the unnormalized inverse transform of forward 
        /// transform.
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x 
        /// contains the transform coeffients.</param>
        public virtual void BackTransform(double[] x)
        {
            Cosqb(ndim, x, wavetable);
        }

        private void Cosqf1(int n, double[] x, double[] wtable)
        {
            int modn, i, k;
            int kc, np2, ns2;
            double xim1;

            ns2 = (n + 1) / 2;
            np2 = n + 2;
            for (k = 1; k < ns2; k++)
            {
                kc = n - k;
                wtable[k + n] = x[k] + x[kc];
                wtable[kc + n] = x[k] - x[kc];
            }
            modn = n % 2;
            if (modn == 0)
            {
                wtable[ns2 + n] = x[ns2] + x[ns2];
            }
            for (k = 1; k < ns2; k++)
            {
                kc = n - k;
                x[k] = wtable[k - 1] * wtable[kc + n] + wtable[kc - 1] * wtable[k + n];
                x[kc] = wtable[k - 1] * wtable[k + n] - wtable[kc - 1] * wtable[kc + n];
            }
            if (modn == 0)
            {
                x[ns2] = wtable[ns2 - 1] * wtable[ns2 + n];
            }
            Rfftf1(n, x, wtable, n);
            for (i = 2; i < n; i += 2)
            {
                xim1 = x[i - 1] - x[i];
                x[i] = x[i - 1] + x[i];
                x[i - 1] = xim1;
            }
        }

        private void Cosqb1(int n, double[] x, double[] wtable)
        {
            int modn, i, k;
            int kc, ns2;
            double xim1;

            ns2 = (n + 1) / 2;
            for (i = 2; i < n; i += 2)
            {
                xim1 = x[i - 1] + x[i];
                x[i] -= x[i - 1];
                x[i - 1] = xim1;
            }
            x[0] += x[0];
            modn = n % 2;
            if (modn == 0)
            {
                x[n - 1] += x[n - 1];
            }
            Rfftb1(n, x, wtable, n);
            for (k = 1; k < ns2; k++)
            {
                kc = n - k;
                wtable[k + n] = wtable[k - 1] * x[kc] + wtable[kc - 1] * x[k];
                wtable[kc + n] = wtable[k - 1] * x[k] - wtable[kc - 1] * x[kc];
            }
            if (modn == 0)
            {
                x[ns2] = wtable[ns2 - 1] * (x[ns2] + x[ns2]);
            }
            for (k = 1; k < ns2; k++)
            {
                kc = n - k;
                x[k] = wtable[k + n] + wtable[kc + n];
                x[kc] = wtable[k + n] - wtable[kc + n];
            }
            x[0] += x[0];
        }

        /// <summary>
        /// Forward cosine FFT with odd wave numbers.
        /// </summary>
        /// <param name="n">Internal parameter.</param>
        /// <param name="x">Internal parameter.</param>
        /// <param name="wtable">Internal parameter.</param>
        protected void Cosqf(int n, double[] x, double[] wtable)
        {
            double sqrt2 = 1.4142135623731;
            double tsqx;

            if (n < 2)
            {
                return;
            }
            if (n == 2)
            {
                tsqx = sqrt2 * x[1];
                x[1] = x[0] - tsqx;
                x[0] += tsqx;
            }
            else
            {
                Cosqf1(n, x, wtable);
            }
        }

        /// <summary>
        /// Backward cosine FFT with odd wave numbers.
        /// </summary>
        /// <param name="n">Internal parameter.</param>
        /// <param name="x">Internal parameter.</param>
        /// <param name="wtable">Internal parameter.</param>
        protected void Cosqb(int n, double[] x, double[] wtable)
        {
            double tsqrt2 = 2.82842712474619;
            double x1;

            if (n < 2)
            {
                x[0] *= 4;
            }
            else if (n == 2)
            {
                x1 = 4 * (x[0] + x[1]);
                x[1] = tsqrt2 * (x[0] - x[1]);
                x[0] = x1;
            }
            else
            {
                Cosqb1(n, x, wtable);
            }
        }

        private void Cosqi(int n, double[] wtable)
        {
            double pih = Math.PI / 2.0D; //1.57079632679491;
            int k;
            double dt;

            dt = pih / (double) n;
            for (k = 0; k < n; k++)
            {
                wtable[k] = Math.Cos((k + 1) * dt);
            }
            Rffti1(n, wtable, n);
        }
    }
}